% Appendix A

\chapter{Program codes} % Main appendix title

\label{AppendixA} % For referencing this appendix elsewhere, use \ref{AppendixA}

\lhead{Appendix A. \emph{Program codes}} % This is for the header on each page - perhaps a shortened title

%-------------------------------------------------------------------
\textbf{SVN Link:} \url{https://code.google.com/p/gpu-em-nim/}

\lstset{language=Matlab, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black},caption=FDTD Free space simulation,label=matlab1}
\begin{lstlisting}
SIZE=201;
ez=zeros(1,SIZE);
hy=zeros(1,SIZE);
imp0=377; %squareroot(u0/e0)
maxTime = 1001;
mm=0;
temp=0;
    for qTime = 1:(maxTime-1)
%        Update Magnetic field
        for  mm = 1:(SIZE - 2)
            hy(mm) = hy(mm) + (ez(mm + 1) - ez(mm)) / imp0;
        end
%         Update Electrical filed
        for mm = 2:(SIZE-1)
            ez(mm) = ez(mm) + (hy(mm) - hy(mm - 1)) * imp0;
        end
%         Source node (hard coded)
        ez(1) = exp(-(qTime - 30) * (qTime - 30) / 100.);
%         Plotting
        figure(1);
        grid on; 
        subplot(2,1,1);
        plot(1:SIZE,ez);
        title('Electirc Component');
        Xlim([0 SIZE]);
        ylim([-1.2 1.2]);
        subplot(2,1,2);
        plot(1:SIZE,hy);
        title('Magnetic Component');
        Xlim([0 SIZE]);
        ylim([-0.003 0.003]);
        pause(0.02);
    end
	
\end{lstlisting}

\lstset{language=Matlab, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black},caption=Perfect magnetic conductor boundary condition,label=matlab2}
\begin{lstlisting}
SIZE=201;
ez=zeros(1,SIZE);
hy=zeros(1,SIZE);
mu=[1.2566e-006*ones(1,SIZE)];   %permeability of free sapce
% epsilon=[8.8542e-012*ones(1,SIZE)]; %permittivity of free space 
epsilon=[8.8542e-012*ones(1,SIZE-80) 1.7708e-011*ones(1,20) 8.8542e-012*ones(1,60)]; %introducing medium of 2*e with a width of 20
% ideal Condition --> Sc= c*delt/delx = 1
delt=1;
delx=299792458;
maxTime = 1001;
mm=0;
temp=0;

    for qTime = 1:(maxTime-1)
%        Update Magnetic field
        for  mm = 1:(SIZE - 2)
            hy(mm) = hy(mm) + (ez(mm + 1) - ez(mm)) * (delt/(delx*mu(mm)));
        end
%         Update Electrical filed
        for mm = 2:(SIZE-1)
            ez(mm) = ez(mm) + (hy(mm) - hy(mm - 1)) * (delt/(delx*epsilon(mm))) ;
        end
%         Source node (hard coded)
        ez(1) = exp(-(qTime - 30) * (qTime - 30) / 100.);
%         Plotting
        figure(1);
        grid on; 
        subplot(2,1,1);
        plot(1:SIZE,ez);
        title('Electirc Component');
        Xlim([0 SIZE]);
        ylim([-1.5 1.5]);
        line([SIZE-80 SIZE-80],[-1.5 1.5],'Color','Red')
        line([SIZE-60 SIZE-60],[-1.5 1.5],'Color','Red')
        subplot(2,1,2);
        plot(1:SIZE,hy);
        title('Magnetic Component');
        Xlim([0 SIZE]);
        ylim([-0.005 0.005]);
        line([SIZE-80 SIZE-80],[-0.005 0.005],'Color','Red')
        line([SIZE-60 SIZE-60],[-0.005 0.005],'Color','Red')
        pause(0.02);
    end
\end{lstlisting}

\lstset{language=Matlab, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black},caption=Absorbing boundary condition (ABC),label=matlab3}
\begin{lstlisting}
clc;
SIZE=201;
ez=zeros(1,SIZE);
hy=zeros(1,SIZE-1);

% Medium Specifications
mu=[1.2566e-006*ones(1,SIZE)];   %permeability of free sapce
% epsilon=[8.8542e-012*ones(1,SIZE)]; %permittivity of free space 
epsilon=[8.8542e-012*ones(1,SIZE-80) 1.7708e-011*ones(1,20) 8.8542e-012*ones(1,60)]; %introducing medium of 2*e with a width of 20

% Courant Number (Accuracy) Sc
% ideal Condition --> Sc= c*delt/delx = 1
delt=1;
delx=299792458;
Sc=299792458*delt/delx;
epsilonr=1;
mur=1;

maxTime = 1001;
mm=0;
temp=0;
ez1q=0;
ez2q=0;
ezmq=0;
ezm1q=0;
    for qTime = 1:(maxTime-1)
%        Update Magnetic field
        for  mm = 1:(SIZE-1)
            hy(mm) = hy(mm) + (ez(mm + 1) - ez(mm)) * (delt/(delx*mu(mm)));
        end
%         Update Electrical filed
        for mm = 2:(SIZE-1)
            ez(mm) = ez(mm) + (hy(mm) - hy(mm - 1)) * (delt/(delx*epsilon(mm))) ;
        end
%         Source node (hard coded)
        ez(2) = ez(2)+exp(-(qTime - 30) * (qTime - 30) / 100.); %additive Source
%         Absorbing Boundary Conditions
        ez(1)=ez2q+(ez(2)-ez1q)*(((Sc/(mur*(epsilonr))^0.5)-1)/((Sc/(mur*(epsilonr))^0.5)+1));
        ez(SIZE)=ezm1q+(ez(SIZE-1)-ezmq)*(((Sc/(mur*(epsilonr))^0.5)-1)/((Sc/(mur*(epsilonr))^0.5)+1));
%         Saving q-1 (pervious step time values) for boundary Conditions
        ez2q=ez(2);
        ez1q=ez(1);
        ezmq=ez(SIZE);
        ezm1q=ez(SIZE-1);
%         Plotting
        figure(1);
        grid on; 
        subplot(2,1,1);
        plot(1:SIZE,ez);
        title('Electirc Component');
        xlim([0 SIZE]);
        ylim([-1.2 1.2]);
        line([SIZE-80 SIZE-80],[-1.2 1.2],'Color','Red') % Medium slab line
        line([SIZE-60 SIZE-60],[-1.2 1.2],'Color','Red') % Medium slab line
        subplot(2,1,2);
        plot(1:SIZE-1,hy);
        title('Magnetic Component');
        xlim([0 SIZE]);
        ylim([-0.005 0.005]);
        line([SIZE-80 SIZE-80],[-0.005 0.005],'Color','Red') % Medium slab line
        line([SIZE-60 SIZE-60],[-0.005 0.005],'Color','Red') % Medium slab line
        pause(0.02);
    end
\end{lstlisting}

\lstset{language=Matlab, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black},caption=Complete 1D FDTD implementation wiht Post Processing,label=matlab4}
\begin{lstlisting}
clc;
SIZE=1024;   
maxTime =1024;

SourceSelect=1; % 0=Sinosoidal, 1=Gauassian

%Constants
c=3e8;

% Courant Number (Accuracy) Sc
% ideal Condition --> Sc= c*delt/delx = 1
% f=3Ghz, lambda=c/f=0.1m, for 4 wavelengths, dx=0.4/(maxtime=1000)
PulseWidth=800;
f=3e9;
w=2*pi*f;    % omega
k0=w/c     ; % free space wave number constant
lambda=c/f;
delx=(4*lambda)/SIZE;
% so dt=dx/c=1.333e-12
delt=delx/c;
Sc=c*delt/delx;
epsilonr=1;
mur=1;


% Incident and Refelected Waves Variables
Eincident=0;
Etransmitted=0;
Etemp=zeros(1,maxTime);


% refractive index variables
Z1 = 750;
z1 = Z1*delx;
Z2 = 760;
z2 = Z2*delx;
Exz1 = zeros(maxTime,1); % record Electric field at 750
Exz2 = zeros(maxTime,1); % record electric field at 760
fspan = 100; % Points to plot in frequency domain

for medium= 1:2
    % Temp Variable
    ez=zeros(1,SIZE);
    hy=zeros(1,SIZE-1);
    mm=0;
    ez1q=0;
    ez2q=0;
    ezmq=0;
    ezm1q=0;
    % Medium Specifications
    mu=1.2566e-006*ones(1,SIZE);   %permeability of free sapce
    if medium==1
        epsilon=8.8542e-012*ones(1,SIZE); % free space
    else
        epsilon=[8.8542e-012*ones(1,SIZE-500) 1.7708e-011*ones(1,500)]; % half medium
    end


    for qTime = 1:(maxTime-1)
%        Update Magnetic field
    	for  mm = 1:(SIZE-1)
            hy(mm) = hy(mm) + (ez(mm + 1) - ez(mm)) * (delt/(delx*mu(mm)));
        end
%         Update Electrical filed
        for mm = 2:(SIZE-1)
            ez(mm) = ez(mm) + (hy(mm) - hy(mm - 1)) * (delt/(delx*epsilon(mm))) ;
        end
        Etemp(qTime)= ez(SIZE-498); %just after boundary of medium
        if SourceSelect==0
%         Source node (hard coded)
		    ez(2) = ez(2)+ (sin(2*pi*(qTime)*f*delt)*Sc);
		else
		    ez(2) = ez(2)+exp(-(qTime - 30) * (qTime - 30) / (PulseWidth./4));
		end
%         Absorbing Boundary Conditions
        ez(1)=ez2q+(ez(2)-ez1q)*(((Sc/(mur*(epsilonr))^0.5)-1)/((Sc/(mur*(epsilonr))^0.5)+1));
        ez(SIZE)=ezm1q+(ez(SIZE-1)-ezmq)*(((Sc/(mur*(epsilonr))^0.5)-1)/((Sc/(mur*(epsilonr))^0.5)+1));
%         Saving q-1 (pervious step time values) for boundary Conditions
		ez2q=ez(2);
		ez1q=ez(1);
		ezmq=ez(SIZE);
		ezm1q=ez(SIZE-1);
%         Plotting
       if medium==2
         if mod(qTime,5)==0   
        figure(1);
        subplot(2,1,1);
        plot(1:SIZE,ez);
        title('Electirc Component');
        xlim([0 SIZE]);
        ylim([-1.2 1.2]);
        if medium==2
            line([SIZE-500 SIZE-500],[-1.2 1.2],'Color','Red') % Medium slab line
        end
        subplot(2,1,2);
        plot(1:SIZE-1,hy);
        title('Magnetic Component');
        xlim([0 SIZE]);
        ylim([-0.005 0.005]);
        if medium==2
        line([SIZE-500 SIZE-500],[-0.005 0.005],'Color','Red') % Medium slab line
        end
         end
       end
          Exz1(qTime)=ez(Z1);
          Exz2(qTime)=ez(Z2);
    end
    if medium==1
        Eincident=Etemp;
    else
        Etransmitted=Etemp;
    end
end



% Frequency Domain Analysis
Fs=1/delt;   %Sampling Frequency
T=1/Fs;      %Sample Time
L=maxTime;      %Length of Signal
time=(0:L-1)*T;    %Time Vector
figure(2);
subplot(2,1,1);
plot(Fs*time,Eincident)
title('Incident Wave');
xlabel('time (picoseconds)')
subplot(2,1,2);
plot(Fs*time,Etransmitted);
title('Transmitted Wave');
xlabel('time (picoseconds)')
% Fourrier Domain
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
FEincident = fft(Eincident,NFFT)/L;
FEtransmitted = fft(Etransmitted,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);          %frequency scaling
% Plot single-sided amplitude spectrum.
figure(3);
subplot(2,1,1);
plot(f,2*abs(FEincident(1:NFFT/2+1))) 
xlim([0 1e11]);
% ylim([0 0.5]);
title('Single-Sided Amplitude Spectrum of Incident Wave')
xlabel('Frequency (Hz)')
ylabel('|Eincident(f)|')
subplot(2,1,2);
plot(f,2*abs(FEtransmitted(1:NFFT/2+1)))
xlim([0 1e11]);
% ylim([0 0.5]);
title('Single-Sided Amplitude Spectrum of Transmitted Wave')
xlabel('Frequency (Hz)')
ylabel('|Etransmitted(f)|')

cTransmitted=FEtransmitted/FEincident
cReflected=1-cTransmitted
eta1=sqrt(1/1);
eta2=sqrt(2/1);
Gamma=(eta2-eta1)/(eta2+eta1)


% eq 33
EXZ1 = fft(Exz1,NFFT)/L;
EXZ2 = fft(Exz2,NFFT)/L;

nFDTD = (1/(1i*k0*(z1-z2))).*log(EXZ2(1:NFFT/2+1)./EXZ1(1:NFFT/2+1));
figure(4);
plot(f(1:fspan), real(nFDTD(1:fspan)));
title('Refractive index re(n)');
xlabel('Frequency (Hz)');
ylabel('re(n)');
line([3e9 3e9],[-15 1.415],'Color','Red')
line([0e9 3e9],[1.415 1.415],'Color','Red')

ReferectiveIndex=(1/(k0*(760-750)*i))*log(FEtransmitted(760)/FEtransmitted(750))

\end{lstlisting}

\lstset{language=Matlab, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black},caption=Drudes Model,label=matlab5}
\begin{lstlisting}

%% Variables Section
clear;
clc;
snapshot=10;
SIZE=1048;   
maxTime =5*SIZE;

SourceSelect=1; % 0=Sinosoidal, 1=Gauassian

%Constants
c=3e8;

% Courant Number (Accuracy) Sc
% ideal Condition --> Sc= c*delt/delx = 1
PulseWidth=800;
f=3e9;      % 3GHz
w=2*pi*f;    % omega
k0=w/c     ; % free space wave number constant
lambda=c/f;
delx=(4*lambda)/SIZE;
% so dt=dx/c=1.333e-12
delt=delx/c;
Sc=c*delt/delx;
epsilonr=-1;
mur=-1;
mu_nought=1.2566e-006;
epsilon_nought=8.8542e-012;
% Drudes model variables
mu_inf=1;
    epsilon_inf=1;
    omega_pe2= 4*pi*pi*f*f*(epsilon_inf-epsilonr);%Plasma frequency electric squared
    omega_pm2= 4*pi*pi*f*f*(mu_inf-mur);%Plasma frequency magnetic squared
    ro_m=0;
    ro_e=0;
% Incident and Refelected Waves Variables
Eincident=0;
Etransmitted=0;
Etemp=zeros(1,maxTime);

% refractive index variables
Z1 = (SIZE/2)+100;
z1 = Z1*delx;
Z2 = (SIZE/2)+110;
z2 = Z2*delx;
Exz1 = zeros(maxTime,1); % record Electric field at 750
Exz2 = zeros(maxTime,1); % record electric field at 760
fspan = 100; % Points to plot in frequency domain
%% Main Loop
for medium= 1:2
    % Medium Specifications
    if medium==1
        mu=1.2566e-006*ones(1,SIZE);   %permeability of free sapce
        epsilon=8.8542e-012*ones(1,SIZE); % free space
        epsilonr=ones(1,SIZE);
        mur=ones(1,SIZE);
    else
        epsilon=[8.8542e-012*ones(1,SIZE-(SIZE/2)) -8.8542e-012*ones(1,(SIZE/4)) 8.8542e-012*ones(1,(SIZE/4))]; % half medium
        mu=[1.2566e-006*ones(1,SIZE-(SIZE/2)) -1.2566e-006*ones(1,(SIZE/4)) 1.2566e-006*ones(1,(SIZE/4))];
        epsilonr=[ones(1,SIZE-(SIZE/2)) -1*ones(1,SIZE/4) ones(1,SIZE/4)];
        mur=[ones(1,SIZE-(SIZE/2)) -1*ones(1,SIZE/4) ones(1,SIZE/4)];
    end
    for  mm = 1:(SIZE)
        omega_pe2= 4*pi*pi*f*f*(epsilon_inf-epsilonr(mm));%Plasma frequency electric squared
        omega_pm2= 4*pi*pi*f*f*(mu_inf-mur(mm));%Plasma frequency magnetic squared
        m_divide=((4*mu_nought*mu_inf)+(mu_nought*omega_pm2*(delt*delt))+(mu_nought*mu_inf*ro_m*(2*delt)));
        am(mm)= 4/m_divide;
        bm(mm)= (ro_m*2*delt)/m_divide;
        cm(mm)= (4*mu_nought*mu_inf)/m_divide;
        dm(mm)= (-mu_nought*omega_pm2*delt*delt)/m_divide;
        em(mm)= (mu_nought*mu_inf*ro_m*2*delt)/m_divide;
        e_divide=((4*epsilon_nought*epsilon_inf)+(epsilon_nought*omega_pe2*(delt*delt))+(epsilon_nought*epsilon_inf*ro_e*(2*delt)));
        ae(mm)= 4/e_divide;
        be(mm)= (ro_e*2*delt)/e_divide;
        ce(mm)= (4*epsilon_nought*epsilon_inf)/e_divide;
        de(mm)= (-epsilon_nought*omega_pe2*delt*delt)/e_divide;
        ee(mm)= (epsilon_nought*epsilon_inf*ro_e*2*delt)/e_divide;
    end
    % Temp Variable
    ez=zeros(1,SIZE);
    ezn_0=zeros(1,SIZE-1);
    ezn_1=zeros(1,SIZE-1);
    hy=zeros(1,SIZE-1);
    hyn_0=zeros(1,SIZE-1);
    hyn_1=zeros(1,SIZE-1);
    by=zeros(1,SIZE);
    byn_0=zeros(1,SIZE);
    byn_1=zeros(1,SIZE);
    dz=zeros(1,SIZE);
    dzn_0=zeros(1,SIZE);
    dzn_1=zeros(1,SIZE);
    mm=0;
    ez1q=0;
    ez2q=0;
    ezmq=0;
    ezm1q=0;
    
    for qTime = 1:(maxTime-1)
        % Drudes model Calculations
        %storing pervious time steps values
        hyn_1=hyn_0;
        hyn_0=hy;
        byn_1=byn_0;
        byn_0=by;
        ezn_1=ezn_0;
        ezn_0=ez;
        dzn_1=dzn_0;
        dzn_0=dz;
%        Update By
        for  mm = 1:(SIZE-1)
            by(mm) = by(mm) + ((ez(mm) - ez(mm+1)) * (delt/(delx)));
        end
%        Update Magnetic field
        for  mm = 1:(SIZE-1)  %changed it from 1 to 2 because of by(0) index problem at 1
            hy(mm) = (am(mm)*(by(mm)-2*byn_0(mm)+byn_1(mm)))+(bm(mm)*(by(mm)-byn_1(mm)))+(cm(mm)*((2*hyn_0(mm))-(hyn_1(mm))))+(dm(mm)*((2*hyn_0(mm))+(hyn_1(mm))))+(em(mm)*(hyn_1(mm)));
        end
%         Update Dz
        for mm = 2:(SIZE-1)
            dz(mm) = dz(mm) + ((hy(mm-1) - hy(mm)) * (delt/delx));
        end
%         Update Electrical filed
        for mm = 2:(SIZE-1)
            ez(mm) = (ae(mm)*(dz(mm)-(2*dzn_0(mm))+dzn_1(mm)))+(be(mm)*(dz(mm)-dzn_1(mm)))+(ce(mm)*((2*ezn_0(mm))-ezn_1(mm)))+(de(mm)*((2*ezn_0(mm))+ezn_1(mm)))+(ee(mm)*ezn_1(mm));
        end
        Etemp(qTime)= ez(SIZE-(SIZE/2)+2); %just after boundary of medium
        if SourceSelect==0
%         Source node (hard coded)
            ez(1) = (sin(2*pi*(qTime)*f*delt)*Sc);
            dz(1)=epsilonr(1)*epsilon_nought*ez(1);
         else
            ez(1)= exp(-(qTime - 30) * (qTime - 30) / (PulseWidth./4));
            dz(1)=epsilonr(1)*epsilon_nought*ez(1);
        end
%         Absorbing Boundary Conditions
        ez(SIZE)=ezm1q+(ez(SIZE-1)-ezmq)*(((Sc/(mur(1)*(epsilonr(1)))^0.5)-1)/((Sc/(mur(1)*(epsilonr(1)))^0.5)+1));
        dz(SIZE)=epsilonr(1)*epsilon_nought; 
%         Saving q-1 (pervious step time values) for boundary Conditions
		ez2q=ez(2);
		ez1q=ez(1);
		ezmq=ez(SIZE);
		ezm1q=ez(SIZE-1);
%         Plotting
       if medium==2
       if mod(qTime,snapshot)==0
        figure(1);
%         subplot(2,1,1);
        plot(1:SIZE,ez);
        title('Electirc Component');
        xlim([0 SIZE]);
%         ylim([-1.2 1.2]);
%         if medium==2
            line([SIZE-(SIZE/2) SIZE-(SIZE/2)],[-2 2],'Color','Red') % Medium slab line
            line([SIZE-(SIZE/4) SIZE-(SIZE/4)],[-2 2],'Color','Red') % Medium slab line
%         end
       end
       end
          Exz1(qTime)=ez(Z1);
          Exz2(qTime)=ez(Z2);
    end
    if medium==1
        Eincident=Etemp;
    else
        Etransmitted=Etemp;
    end
end
%% Post processing
% Frequency Domain Analysis
Fs=1/delt;   %Sampling Frequency
T=1/Fs;      %Sample Time
L=maxTime;      %Length of Signal
time=(0:L-1)*T;    %Time Vector
figure(2);
subplot(2,1,1);
plot(Fs*time,Eincident)
title('Incident Wave');
xlabel('time (picoseconds)')
subplot(2,1,2);
plot(Fs*time,Etransmitted);
title('Transmitted Wave');
xlabel('time (picoseconds)')
% Fourrier Domain
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
FEincident = fft(Eincident,NFFT)/L;
FEtransmitted = fft(Etransmitted,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);          %frequency scaling
% Plot single-sided amplitude spectrum.
figure(3);
subplot(2,1,1);
plot(f,2*abs(FEincident(1:NFFT/2+1))) 
xlim([0 1e11]);
% ylim([0 0.5]);
title('Single-Sided Amplitude Spectrum of Incident Wave')
xlabel('Frequency (Hz)')
ylabel('|Eincident(f)|')
subplot(2,1,2);
plot(f,2*abs(FEtransmitted(1:NFFT/2+1)))
xlim([0 1e11]);
% ylim([0 0.5]);
title('Single-Sided Amplitude Spectrum of Transmitted Wave')
xlabel('Frequency (Hz)')
ylabel('|Etransmitted(f)|')

cTransmitted=FEtransmitted/FEincident
cReflected=1-cTransmitted
eta1=sqrt(1/1);
eta2=sqrt(-1/1);
Gamma=(eta2-eta1)/(eta2+eta1)


% eq 33
EXZ1 = fft(Exz1,NFFT)/L;
EXZ2 = fft(Exz2,NFFT)/L;

nFDTD = (1/(1i*k0*(z1-z2))).*log(EXZ2(1:NFFT/2+1)./EXZ1(1:NFFT/2+1));
figure(4);
%plot(f(1:fspan), real(nFDTD(1:fspan)));
plot(f(20:50), real(nFDTD(20:50)));
title('Refractive index re(n)');
xlabel('Frequency (Hz)');
ylabel('re(n)');
line([3e9 3e9],[-4 -1],'Color','Red')
line([1e9 3e9],[-1 -1],'Color','Red')

ReferectiveIndex=(1/(k0*(Z2-Z1)*i))*log(FEtransmitted(Z2)/FEtransmitted(Z1))

f0=3e9;
w0=2*pi*f0;
einf=1;
wpe2= 2*(w0^2);
f = 2e9:1e8:4e9;
length(f)

er = einf - wpe2./(2*pi*f).^2;
figure(4)
hold on
plot(f,er)
line([3e9 3e9],[-15 -1],'Color','Red')
line([0e9 3e9],[-1 -1],'Color','Red')

\end{lstlisting}




\lstset{language=[ISO]C++, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black},caption=c++ Implementation with Dynamic file handling, label=cplus1}
\begin{lstlisting}
#include <iostream>
#include <fstream>
#include <sstream>
#include <windows.h>

using namespace std;
const double pi = 3.14;

void main()
{
	// -------- Save to file Variables -------- 
	fstream snapshot;
	std::string filename ;
	std::stringstream stream;
	stream<<"results";
	CreateDirectory(stream.str().c_str(), NULL) ;		//create directory of results

	int SIZE=1000;
    int maxTime = 1024;
	int SourceSelect = 1; 			// 0=Sinosoidal, 1=Gauassian
	cout<<"----Select Source----"<<endl;
	cout<<"0)Sinosoidal 1) Gauassian"<<endl;
	do{
		cin>>SourceSelect;
	}while(SourceSelect != 1 && SourceSelect != 0);
	if (SourceSelect == 0)
		maxTime = 5001;

	// -------- Constants -------- 
	double c = 3e8;

 	// Courant Number (Accuracy) Sc
	// ideal Condition --> Sc= c*delt/delx = 1
	// f=3Ghz, lambda=c/f=0.1m, for 4 wavelengths, dx=0.4/(maxtime=1000)
	double PulseWidth = 800;
	double f = 3e9;
	double w = 2 * pi * f;    		// omega
	double k0 = w/c     ; 			// free space wave number constant
	double lambda = c / f;
	double delx = (4*lambda) / SIZE;
	double delt = delx / c;
	double Sc = c * delt / delx;
	int epsilonr = 1;
	int mur = 1;

    // -------- Incident and Refelected Waves Variables -------- 
	double *Eincident = new double[maxTime] ; 
	for(int j=0; j<maxTime; j++)
		Eincident[j] = 0;
	
	double *Etransmitted = new double [maxTime] ;
	for(int j=0; j<maxTime; j++)
		Etransmitted[j] = 0;

	double *Etemp = new double [maxTime];
	for(int j=0; j<maxTime; j++)
		Etemp[j] = 0;

	// -------- refractive index variables -------- 
	int Z1 = 750;
	double z1 = Z1*delx;
	int Z2 = 760;
	double z2 = Z2*delx;

	double *Exz1 =new double [maxTime];
	for (int i=0; i<maxTime ; i++)
		Exz1[i] = 0;

	double *Exz2 =new double [maxTime];
	for (int i=0; i<maxTime ; i++)
		Exz2[i] = 0;
	
	double *mu = new double [SIZE];
	for(int j=0; j<SIZE; j++)
			mu[j] = 1.2566e-006;

	double *epsilon = new double [SIZE];

	double *ez = new double [SIZE];
	for(int j=0; j<SIZE; j++)
		ez[j] = 0;

	double *hy = new double [SIZE-1];
	for(int j=0; j<SIZE-1; j++)
		hy[j] = 0;
	
	for (int medium=1; medium<=2; medium++)
	{
    	// Temp Variable
		int mm = 0;
		double ez1q = 0;
		double ez2q = 0;
		double ezmq = 0;
		double ezm1q = 0;
		
		// -------- Medium Specifications -------- 
		if (medium==1)
		{
			cout<<"Calculating Wave propagation in Free Space"<<endl;
			for(int j=0; j<SIZE; j++)
				epsilon[j] = 8.8542e-012;   						//epsilon=8.8542e-012*ones(1,SIZE); // free space
		}
		else
		{
			cout<<"Calculating Wave propagation in denser Medium"<<endl;
			int j;
			for(j=0; j<SIZE-500; j++)								//epsilon=[8.8542e-012*ones(1,SIZE-500) 1.7708e-011*ones(1,500)]; // half medium
				epsilon[j] = 8.8542e-012;
			for(int k=j ;k<SIZE ;k++)
				epsilon[j] = 1.7708e-011;
		}

		// -------- Main Loop -------- 
		for (int qTime=0; qTime<(maxTime-1); qTime++)
		{
			for(int nn=0; nn<(SIZE-1); nn++)
			{
				hy[nn] = hy[nn] + (ez[nn+1] - ez[nn]) * (delt/(delx*mu[nn]));				//Update Magnetic field
			}

			for(int nn=1; nn<(SIZE-1); nn++)
			{
				ez[nn] = ez[nn] + (hy[nn] - hy[nn-1]) * (delt/(delx*epsilon[nn]));			//Update Electrical field
			}

			if (SourceSelect==0)															//Source node
			    ez[1] = ez[1] + (sin(2*pi*(qTime)*f*delt)*Sc);
			else
			    ez[1] = ez[1] + exp((-(qTime+1 - 30) * (qTime - 30)) / (PulseWidth/4));

	        Etemp[qTime]= ez[SIZE-498]; 													//Save ez after boundary
		    
		    // -------- Absorbing Boundary Conditions -------- 
	        ez[0] = ez2q+(ez[1]-ez1q)*( ((Sc/pow(mur*epsilonr,0.5))-1 ) / ((Sc/pow(mur*epsilonr,0.5))+1) );	
	        ez[SIZE-1] = ezm1q+(ez[SIZE-1-1]-ezmq)*( ((Sc/pow(mur*epsilonr,0.5))-1 ) / ((Sc/pow(mur*epsilonr,0.5))+1) );
		    
		    // -------- Saving pervious step time values -------- 
			ez2q = ez[1]; 
			ez1q = ez[0]; 
			ezmq = ez[SIZE-1];
			ezm1q= ez[SIZE-2];
			Exz1[qTime] = ez[Z1-1];
	        Exz2[qTime] = ez[Z2-1];

	        // -------- Saving to file -------- 
			stream.str(std::string());   						// clear stringstream
			stream<<"./results/"<<"Efield"<<medium<<"_"<<qTime<<".jd";   		// concatenate
			filename = stream.str();		 					// copy string
			snapshot.open(filename, ios::out|ios::binary);
			for (mm = 0; mm < SIZE; mm++)
				snapshot.write((char *)&ez[mm],sizeof(double));
			snapshot.close();
		} //end qTime

		if (medium==1)
		{
			for(int i=0;i<maxTime;i++)
			Eincident[i] = Etemp[i];
		}
		else
		{
			for(int i=0;i<maxTime;i++)
			Etransmitted[i] = Etemp[i];
		}
		
	} //end medium loop
	
	cout<<"Writing Values to files"<<endl;
	stream.str(std::string());
	stream<<"./results/"<<"Eincident"<<".jd";
	filename = stream.str();
	snapshot.open(filename, ios::out|ios::binary);
	for (int mm = 0; mm < SIZE; mm++)
	snapshot.write((char *)&Eincident[mm],(sizeof(double)));
	snapshot.close();

	stream.str(std::string());
	stream<<"./results/"<<"Etransmitted"<<".jd";
	filename = stream.str();
	snapshot.open(filename, ios::out|ios::binary);
	for (int mm = 0; mm < SIZE; mm++)
	snapshot.write((char *)&Etransmitted[mm],(sizeof(double)));
	snapshot.close();

	stream.str(std::string());
	stream<<"./results/"<<"Exz1"<<".jd";
	filename = stream.str();
	snapshot.open(filename, ios::out|ios::binary);
	for (int mm = 0; mm < SIZE; mm++)
	snapshot.write((char *)&Exz1[mm],(sizeof(double)));
	snapshot.close();

	stream.str(std::string());
	stream<<"./results/"<<"Exz2"<<".jd";
	filename = stream.str();
	snapshot.open(filename, ios::out|ios::binary);
	for (int mm = 0; mm < SIZE; mm++)
	snapshot.write((char *)&Exz2[mm],(sizeof(double)));
	snapshot.close();

	stream.str(std::string());
	stream<<"./results/"<<"maxTime"<<".jd";
	filename = stream.str();
	snapshot.open(filename, ios::out|ios::binary);
	snapshot.write((char *)&maxTime,(sizeof(int)));
	snapshot.close();

	stream.str(std::string());
	stream<<"./results/"<<"data"<<".jd";
	filename = stream.str();
	snapshot.open(filename, ios::out|ios::binary);
	snapshot.write((char *)&delt,(sizeof(double)));
	snapshot.write((char *)&k0,(sizeof(double)));
	snapshot.write((char *)&z1,(sizeof(double)));
	snapshot.write((char *)&z2,(sizeof(double)));
	snapshot.close();

	cout<<endl;
	cout<<"Enter anything to exit"<<endl;
	cin>>SourceSelect;

	// -------- Memory Deallocaion -----------
	delete[] Eincident;
	delete[] Etransmitted;
	delete[] Etemp;
	delete[] mu;
	delete[] epsilon;
	delete[] ez;
	delete[] hy;
	delete[] Exz1;
	delete[] Exz2;
}
\end{lstlisting}

\lstset{language=Matlab, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black},caption=FDTD Free space simulation,label=cplus2}
\begin{lstlisting}
SIZE=1001;
medium=1;
fspan=100;
filepath=fullfile(pwd, 'results');

fidp = fopen(strcat(filepath,'\maxTime','.jd'),'r','l');
if fidp==-1
    display('Error');
    return
end
maxTime=fread(fidp,1,'int');
fclose(fidp);

fidp = fopen(strcat(filepath,'\data','.jd'),'r','l');
if fidp==-1
    display('Error');
    return
end
data=fread(fidp,4,'double');
fclose(fidp);
delt=data(1);
k0=data(2);
z1=data(3);
z2=data(4);

fidp = fopen(strcat(filepath,'\Eincident','.jd'),'r','l');
if fidp==-1
    display('Error');
    return
end
Eincident=fread(fidp,maxTime,'double');
fclose(fidp);

fidp = fopen(strcat(filepath,'\Etransmitted','.jd'),'r','l');
if fidp==-1
    display('Error');
    return
end
Etransmitted=fread(fidp,maxTime,'double');
fclose(fidp);

% Frequency Domain Analysis
Fs=1/delt;   %Sampling Frequency
T=1/Fs;      %Sample Time
L=maxTime;      %Length of Signal
time=(0:L-1)*T;    %Time Vector
figure(1);
subplot(2,1,1);
plot(Fs*time,Eincident)
title('Incident Wave');
xlabel('time (picoseconds)')
subplot(2,1,2);
plot(Fs*time,Etransmitted);
title('Transmitted Wave');
xlabel('time (picoseconds)')
% Fourrier Domain
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
FEincident = fft(Eincident,NFFT)/L;
FEtransmitted = fft(Etransmitted,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);          %frequency scaling
% Plot single-sided amplitude spectrum.
figure(2);
subplot(2,1,1);
plot(f,2*abs(FEincident(1:NFFT/2+1))) 
xlim([0 1e11]);
% ylim([0 0.5]);
title('Single-Sided Amplitude Spectrum of Incident Wave')
xlabel('Frequency (Hz)')
ylabel('|Eincident(f)|')
subplot(2,1,2);
plot(f,2*abs(FEtransmitted(1:NFFT/2+1)))
xlim([0 1e11]);
% ylim([0 0.5]);
title('Single-Sided Amplitude Spectrum of Transmitted Wave')
xlabel('Frequency (Hz)')
ylabel('|Etransmitted(f)|')
 
cTransmitted=FEtransmitted/FEincident;
cReflected=1-cTransmitted;
eta1=sqrt(1/1);
eta2=sqrt(2/1);
Gamma=(eta2-eta1)/(eta2+eta1) 
 
% eq 33
fidp = fopen(strcat(filepath,'\Exz1','.jd'),'r','l');
if fidp==-1
    display('Error');
    return
end
Exz1=fread(fidp,maxTime,'double');
fclose(fidp);

fidp = fopen(strcat(filepath,'\Exz2','.jd'),'r','l');
if fidp==-1
    display('Error');
    return
end
Exz2=fread(fidp,maxTime,'double');
fclose(fidp);

EXZ1 = fft(Exz1,NFFT)/L;
EXZ2 = fft(Exz2,NFFT)/L;

nFDTD = (1/(1i*k0*(z1-z2))).*log(EXZ2(1:NFFT/2+1)./EXZ1(1:NFFT/2+1));
figure(3);
plot(f(1:fspan), real(nFDTD(1:fspan)));
title('Refractive index re(n)');
xlabel('Frequency (Hz)');
ylabel('re(n)');
line([3e9 3e9],[-15 1.415],'Color','Red')
line([0e9 3e9],[1.415 1.415],'Color','Red')

for medium= 1:2
    j=0;
    for qTime = 1:(maxTime-1)
        fidp = fopen(strcat(filepath,'\Efield',int2str(medium),'_',int2str(j),'.jd'),'r','l');
        if fidp==-1
            display('Error');
            return
        end
        ez=fread(fidp,maxTime,'double');
        fclose(fidp);
        j=j+1;
% Plotting
        figure(4);
        subplot(2,1,1);
        plot(1:SIZE,ez);
        title('Electirc Component');
        xlim([0 SIZE]);
        ylim([-1.2 1.2]);
        if medium==2
            line([SIZE-500 SIZE-500],[-1.2 1.2],'Color','Red') % Medium slab line
        end
        subplot(2,1,2);
%         plot(1:SIZE-1,hy);
        title('Magnetic Component');
        xlim([0 SIZE]);
        ylim([-0.005 0.005]);
        if medium==2
        line([SIZE-500 SIZE-500],[-0.005 0.005],'Color','Red') % Medium slab line
        end
    end
end
\end{lstlisting}

\lstset{language=[ISO]C++, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black},caption=Kernel function in OpenCL, label=gpu22}
\begin{lstlisting}

#define PRECISION float

/*
#ifdef cl_amd_fp64
#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#else
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#endif
*/

__kernel void hy_kernel(__global PRECISION *hy, __global PRECISION *ez, __global PRECISION *mu, const PRECISION delt, const PRECISION delx, const int SIZE) 
{
	unsigned int i = get_global_id(0);
  if(i < SIZE)
    hy[i] = hy[i] + (ez[i+1] - ez[i]) * (delt/(delx*mu[i]));
}
__kernel void ez_kernel(__global PRECISION *hy, __global PRECISION *ez, __global PRECISION *epsilon, const PRECISION delt, const PRECISION delx, const int SIZE,__global PRECISION *Etemp,__global PRECISION *Exz1,__global PRECISION *Exz2,const int qTime,const PRECISION mur,const PRECISION epsilonr,const PRECISION Sc,PRECISION ez1q,PRECISION ez2q,PRECISION ezmq,PRECISION ezm1q,const PRECISION pi,const PRECISION f,const int SourceSelect,const PRECISION PulseWidth,const int Z1,const int Z2)
{
	unsigned int i = get_global_id(0);
  if(i > 0 && i < SIZE)
    ez[i] = ez[i] + (hy[i] - hy[i-1]) * (delt/(delx*epsilon[i]));

  if (SourceSelect==0)	
  {														//Source node
	ez[1] = ez[1] + (sin(2*pi*(qTime)*f*delt)*Sc);
	}
  else
  {
	ez[1] = ez[1] + exp((-(qTime+1 - 30) * (qTime - 30)) / (PulseWidth/4));
	}
  Etemp[qTime]= ez[SIZE-(SIZE/2)+2]; 													//Save ez after boundary
  // -------- Absorbing Boundary Conditions -------- 
  ez[0] = ez2q+(ez[1]-ez1q)*( ((Sc/pow(mur*epsilonr,0.5f))-1 ) / ((Sc/pow(mur*epsilonr,0.5f))+1.f) );
  ez[SIZE-1] = ezm1q+(ez[SIZE-1-1]-ezmq)*( ((Sc/pow(mur*epsilonr,0.5f))-1.f ) / ((Sc/pow(mur*epsilonr,0.5f))+1.f) );	  
  
  barrier(CLK_GLOBAL_MEM_FENCE|CLK_LOCAL_MEM_FENCE);    
  // -------- Saving pervious step time values -------- 
  ez2q = ez[1]; 
  ez1q = ez[0]; 
  ezmq = ez[SIZE-1];
  ezm1q= ez[SIZE-2];
  Exz1[qTime] = ez[Z1-1];
  Exz2[qTime] = ez[Z2-1];
}

\end{lstlisting}